
!-------------------------------------------------------------------------------------
subroutine evolve_bgk(Fbgk, Ebgk, dt)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit
  use fft_wrap

  implicit none

  real,    intent(inout) :: Fbgk(jv, jz)
  real,    intent(out)   :: Ebgk(jz)
  real,    intent(in)    :: dt

  integer                :: iz, iv

  complex                :: fz(jz)
  complex                :: fv(jv)


  !-- advect in z-direction

  do iv=1,jv
     fz(:) = Fbgk(iv,:)
     call fft_forward(jz,fz,fz)
     fz = fz/float(jz)
     fz(:) = fz(:) * exp ( -iunit * kZvec(:) * dt * vVec(iv) )
     call fft_backward(jz,fz,fz)
     Fbgk(iv,:) = fz(:)
  enddo

  call efield_bgk(Fbgk, Ebgk)


  !-- advect in v-direction

   do iz=1,jz
     fv(:) = Fbgk(:,iz)
     call fft_forward(jv,fv,fv)
     fv = fv/float(jv)
     fv(:) = fv(:) * exp ( -iunit * kVvec(:) * dt * Ebgk(iz) )
     call fft_backward(jv,fv,fv)
     Fbgk(:,iz) = fv(:)
  enddo


end subroutine evolve_bgk


!-------------------------------------------------------------------------------------
subroutine efield_bgk(Fbgk, Ebgk)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit
  use fft_wrap

  implicit none

  real,    intent(in)    :: Fbgk(jv, jz)
  real,    intent(inout) :: Ebgk(jz)

  integer                :: iz, iv

  complex                :: fz(jz)
  complex                :: fv(jv)


  !-- compute charge density; store in fz

  do iz=1,jz  
     fz(iz) = sum(Fbgk(:,iz))*dv 
  end do


  !-- find electric field

  call fft_forward(jz,fz,fz)
  fz = fz/float(jz)
  fz(:) = - iunit*fz(:)/kZvec(:)
  fz(1) = 0;
  call fft_backward(jz,fz,fz)
  Ebgk  = fz;


end subroutine efield_bgk


!-------------------------------------------------------------------------------------
subroutine derivV_bgk(Fbgk, DVbgk)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit

  use fft_wrap


  implicit none

  real,    intent(in)    :: Fbgk(jv, jz)
  real,    intent(out)   :: DVbgk(jv, jz)

  integer                :: iz
  complex                :: fv(jv)


    !-- take derivative in v-direction

    do iz=1,jz
       fv(:) = Fbgk(:,iz)
       call fft_forward(jv,fv,fv)
       fv = fv/float(jv)
       fv(:) = iunit*fv(:)*kVvec(:)
       call fft_backward(jv,fv,fv)
       DVbgk(:,iz) = fv(:)
    enddo


end subroutine derivV_bgk




!-------------------------------------------------------------------------------------
subroutine derivZ_bgk(Fbgk, DZbgk)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit

  use fft_wrap


  implicit none

  real,    intent(in)    :: Fbgk(jv, jz)
  real,    intent(out)   :: DZbgk(jv, jz)

  integer                :: iv
  complex                :: fz(jz)



    !-- take derivative in z-direction

    do iv=1,jv
       fz(:) = Fbgk(iv,:)
       call fft_forward(jz,fz,fz)
       fz = fz/float(jz)
       fz(:) = iunit*fz(:)*kZvec(:)
       call fft_backward(jz,fz,fz)
       DZbgk(iv,:) = fz(:)
    enddo


end subroutine derivZ_bgk



!-------------------------------------------------------------------------------------
subroutine sumUz_bgk(Fbgk,  Ubgk)
!-------------------------------------------------------------------------------------

  use dimensions
  implicit none

  real,    intent(in)    :: Fbgk(jv, jz)
  real,    intent(out)   :: Ubgk(jz)

  integer    :: iz


  !-- integrate over "v" to find longtitudinal velocity Uz -- 
       
  do iz = 1,jz 
      Ubgk(iz) = sum( Fbgk(:,iz) * vVec ) / sum( Fbgk(:,iz) )
  enddo
   

end subroutine sumUz_bgk



!-------------------------------------------------------------------------------------
subroutine bgk_test_equil(Fbgk, DVbgk, Ebgk)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: runname, dt, nt

  use mpi_wrap


  implicit none

  real,    intent(inout) :: Fbgk(jv, jz)
  real,    intent(out)   :: DVbgk(jv, jz)
  real,    intent(out)   :: Ebgk(jz)

  real                   :: DZbgk(jv, jz)
  real                   :: equil1(jv, jz)
  real                   :: equil2(jv, jz)
  real                   :: rhs(jv, jz)

  integer                :: iv, iz, it



  !-- evolve toward equilibrium

  it = 0
 
  do while (it .le. nt) 

    call evolve_bgk(Fbgk, Ebgk, dt)

    call efield_bgk(Fbgk, Ebgk)

    call derivV_bgk(Fbgk, DVbgk)

    call derivZ_bgk(Fbgk, DZbgk)

    !-- check if distribution is in equilibrium

    do iv=1,jv
       equil1(iv,:) = - vvec(iv)*DZbgk(iv,:)
    end do
 
    do iz=1,jz
       equil2(:,iz) = - Ebgk(iz)*DVbgk(:,iz)
    end do

    rhs = equil1 + equil2

    !Fbgk = Fbgk + rhs*dt
 
    it = it + 1


    write(*,*) it*dt, maxval(abs(Fbgk)), maxval(abs(rhs)), maxval(abs(rhs/Fbgk))

end do



  open(91, file= trim(runname)//'.df', form='unformatted', access='stream')
  write(91) rhs
  close(91)

  open(92, file= trim(runname)//'.f', form='unformatted', access='stream')
  write(92) Fbgk
  close(92)


!  call mpi_allreduce(tmp, dens, jz*jy*jx, &
!      MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, mpi_err)


end subroutine bgk_test_equil

!-------------------------------------------------------------------------------------


